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We investigate the Hubbard model on a two-dimensional square lattice by the perturbation 
expansion to the fourth order in the on-site Coulomb repulsion U. Numerically calculating 
all diagrams up to the fourth order in self-energy, we examine the convergence of perturbation 
series in the lattice system. We indicate that the coefficient of each order term rapidly decreases 
as in the impurity Anderson model for T > O.li in the half-filled case, but it holds in the doped 
case even at lower temperatures. Thus, we can expect that the convergence of perturbation 
expansion in U is very good in a wide parameter region also in the lattice system, except 
for T < O.lt in the half-filled case. We next calculate the density of states in the fourth-order 
perturbation. In the half-filled case, the shape in a moderate correlation regime is quite different 
from the three peak structure in the second-order perturbation. Remarkable upper and lower 
Hubbard bands locate at uj ~ ±?7/2, and a pseudogap appears at the Fermi level ui = 0. This 
is considered as the precursor of the Mott-Hubbard antiferromagnetic structure. In the doped 
case, quasiparticles with very heavy mass are formed at the Fermi level. Thus, we conclude 
that the fourth-order perturbation theory overall well explain the asymptotic behaviors in a 
strong correlation regime. 
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1. Introduction 

Over the past several decades, research on high-T c 
cuprates has been a major theme in condensed mat- 
ter physics. Unconventional superconductivity with line 
nodes, not an isotropic gap in the BCS theory, appears 
with carrier doping in the Mott insulator. Inevitably, a 
strong correlation has been considered very important. 
The ground state in the Mott insulator is the antifer- 
romagnetic (AF) state, and the interplay between mag- 
netism and unconventional superconductivity has been 
the key issue. In order to understand the phase diagram 
and physical properties, a large number of theoretical and 
experimental research studies have been performed. 1 ' 2 
Today, we understand that superconductivity originates 
from the AF spin fluctuation. On the same footing, most 
anomalous features in cuprates can be consistently ex- 
plained. 3,4 

Quantitatively, fluctuation-exchange approximation 
(FLEX) has been used as an efficient numerical method. 
This method successfully describes the behavior of the 
Fermi liquid (FL) state in the vicinity of the AF criti- 
cal point on the basis of the microscopic Hamiltonian. 
However, it has a few disadvantages. The specific mode 
of the spin fluctuation is overestimated owing to the 
enhancement by the partial summation of bubble and 
ladder diagrams. The Hubbard peak structure in the 
density of states (DOS) is smeared and unclear. The 
Mott transition cannot be described. These failures in 
the normal state are concerned with the fact that the 
Mott-Hubbard character arises from a local correlation. 
The essence has been clarified by dynamical mean-field 



theory (DMFT). 5,6 Nowadays, great efforts have been 
made to take short-range correlations into account, 7-13 
such as the dynamical cluster approximation (DC A). 
As for superconductivity, on the other hand, the FLEX 
well explains the d-wave spin-singlet superconductivity 
near the AF phase, and the p-wave spin-triplet channel 
can be dominant near the ferromagnetic phase. These 
are reasonable results. However, the FLEX cannot sim- 
ply explain the p-wave spin-triplet superconductivity in 
Sr2Ru04. In order to obtain the triplet state in a single- 
band Hubbard model for the 7 band, which is consid- 
ered to trigger the superconductivity in Sr2Ru04, an ad- 
ditional off-site interaction is necessary at least. 14 This 
may be related to the fact that the spin fluctuation in 
the 7 band shows a featureless broad hump at around 
the r-ma point, although the FLEX is justified in sys- 
tems with remarkable spin fluctuation. There is a pos- 
sibility that the FLEX ignores crucial terms in the p- 
wave state in Sr 2 Ru04. In fact, the third-order per- 
turbation theory 15, 16 and the one-loop renormalization 
group method 17, 18 indicate a p-wave state within a con- 
ventional single-band Hubbard model for the 7 band. 
In these cases, the important process in the pairing in- 
teraction mainly originates from the third-order term. 
This is the case for the spin-triplet channel in the two- 
dimensional electron gas system. 19,20 In the FLEX, the 
process, which is classified into a vertex correction term, 
is not included. Thus, the FLEX may not work except in 
the vicinity of the critical point such as the AF (or ferro- 
magnetic) phase, although we need further investigation 
to clarify it. 

On the other hand, the third-order perturbation the- 
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ory seems to be efficient in the study of superconductiv- 
ity. It is comprehensively applicable in both spin singlet 
and triplet states. Input parameters are only dispersion 
relation and carrier number. However, this method is an 
approach based on a weak correlation. The applicabil- 
ity in strongly correlated systems is obscure. We need 
to evaluate the convergence of the perturbation expan- 
sion. In addition, we can expect that a very large mass 
enhancement factor is obtained in higher-order perturba- 
tion from previous studies. 21, 22 Thus, higher-order per- 
turbation in the on-site Coulomb repulsion U is impor- 
tant not only as fundamental knowledge in strongly cor- 
related systems, but also as a quantitative method. Up to 
now, from numerical calculations restricted to k points 
on the Fermi surface (FS), Nomura and Yamada have 
discussed the convergence of the pairing interaction in 
the fourth-order perturbation, 23 and Shinkai et. al have 
evaluated the mass renormalization factor. 21 ' 22 However, 
owing to the restriction, physical properties in the fourth- 
order perturbation have not been clarified yet. Thus, in 
this study, we carry out the perturbation expansion to 
the fourth order in the entire first Brillouin zone. Its cal- 
culation is very instructive, and clarifies physical proper- 
ties that are indefinite in previous research studies. 21-24 

In the present study, one of theoretical backgrounds 
is the result of the perturbation approach of Yamada 
and Yosida for the impurity Anderson model. 25-27 They 
carried out the perturbation expansion to the fourth or- 
der in U in the context of the Kondo problem. Today, 
we have the exact solution, and the applicability of the 
perturbation theory has been confirmed. 28 The radius of 
convergence is infinite, and physical quantities are ana- 
lytic in U. For instance, the spin susceptibility and the 
specific heat coefficient are enhanced like the exponen- 
tial function as a function of U. Thus, the coefficient of 
each nth-order term in these physical quantities rapidly 
decreases almost in proportion to ~ 1/n!. Because of 
this remarkable property, even if we truncate the per- 
turbation expansion at a finite order, physical quantities 
rapidly approach the exact values with increasing cutoff 
order such as 2, 4, 6, • • • . For instance, the exact Wil- 
son ratio is ~ 1.962 at u = U/irAta = 2. 28 This can 
be regarded as a sufficiently strong correlation regime, 
since the exact value is 2 at a strong correlation limit. 
In this case, the approximate values are ~ 1.639 for the 
second order, ~ 1.889 for the fourth order and ~ 1.952 
for the sixth order. This indicates that the fourth-order 
perturbation expansion has sufficient accuracy in a mod- 
erate correlation regime. Such good convergence is one 
of characteristics in the FL state. 

Since we have no exact solution in two- or higher- 
dimensional lattice systems, we cannot guarantee the 
convergence of the perturbation expansion. In fact, the 
ground states in many lattice systems are not the FL 
state, but the magnetic or superconducting state. We 
will need much higher-order perturbation terms to re- 
store the critical fluctuation near these critical points. 
It may be rather better to perform partial summations, 
such as the FLEX. However, above these transition tem- 
peratures, the system can be considered to be in the FL 
state, because of the principle of adiabatic continuation 
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stressed by Anderson. 29 As long as no phase transition 
occurs, the system connects adiabatically with the non- 
interacting system, and stays in the FL state. In this 
case, the perturbation expansion is still applicable. The 
physical quantities of the system will asymptotically ap- 
proach the exact behavior with increasing cutoff order 
in the perturbation expansion. In this paper, we perform 
the perturbation expansion for the normal self-energy. 
Since the first-order term provides only a constant shift, 
the second-order term is the first one that includes cor- 
relation effects due to the on-site Coulomb repulsion U. 
The next significant term comes from the fourth-order 
term. This is because the third-order term is relatively 
small owing to the fact that it vanishes in systems with 
particle-hole symmetry. The fourth-order term can in- 
clude correlation effects that have not been grasped in 
the third-order perturbation theory so far. This can qual- 
itatively change the asymptotic behavior in the strong 
correlation regime. In addition, we can estimate the va- 
lidity of the perturbation expansion to the third order 
by comparing each order term. If the perturbation ex- 
pansion has good convergence of 1/n! as in the impurity 
case, the fourth-order perturbation theory will be valid 
in a wider parameter region than the third-order per- 
turbation. Thus, the fourth-order perturbation theory in 
lattice systems can be considered as one of several ef- 
ficient methods of studying strongly correlated systems. 
The investigation of the fourth-order perturbation is very 
fruitful not only as fundamental knowledge in strongly 
correlated systems, but also as a quantitative method. 

In this paper, we calculate the fourth-order term for 
the normal self-energy in lattice systems, and examine 
the convergence of the perturbation expansion. In addi- 
tion, we investigate single-particle quantities, such as the 
DOS, in the fourth-order perturbation theory. In particu- 
lar, we clarify the asymptotic behavior in the strong cor- 
relation regime. It is meaningful that we compare these 
results with those given by many other theoretical stud- 
ies. Furthermore, it will be useful to introduce numerical 
algorithms to evaluate vertex correction terms, which do 
not have the convolution form. 

This paper is organized as follows. In §2, we briefly 
formalize the perturbation expansion to the fourth order 
for the self-energy. In §3, we describe some techniques 
in numerical calculations, particularly how to estimate 
vertex correction terms, which are usually neglected be- 
cause of not being in the convolution form. We present 
numerical results in both the half-filled and doped cases 
in §4. We discuss the convergence of the perturbation 
expansion in the former half, and then show properties 
of single-particle quantities in the latter half. Finally, in 
§5, we summarize our study and give an outline of future 
works. 

2. Formalism 

In this paper, we study the Hubbard Hamiltonian on 
an N x N square lattice: 

ka i 
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U is the on-site Coulomb repulsion, and \i is the chemical 
potential. The dispersion relation efc is given by 

e fc = -2i (cos(fcr) + cos(ky)) , (2) 

where the wave vector k = (k x ,k y ) is measured in units 
of the inverse of lattice constant. The single-particle 
Green's function is expressed with the self-energy S(k) 
as 

g(k)- 1 = go(k)- 1 - s(k) = iw n - & - s(k), (3) 

where = — ^ and k = (k,iu) n ). w n = (2n + l)irT 
is the fermion Matsubara frequency at a temperature T. 
Here, we set the Boltzmann constant k m athrmB = 1. 
Hereafter, we measure energy in units of t, that is, wc 
set t = 1. For a given chemical potential /j,, the electron 
density n is determined by 

» = £ m = £ s°( fc ) = ^2 E w 

<7 /ccr /co- 

in the noninteracting system, and 

n = e aw = E (^w - ^( fc )) + ^2 E /(&). (5) 

in the interacting system. In these final expressions, 

^) = e WrT = K 1 - tanh ^))' (6) 

is the Fermi-Dirac distribution function with (3 = 1/T, 
and the sum over k denotes 

E — ~jy2 E E ' 

k k iuj n 

The final expressions in eqs. (5) and (6) are convenient 
in numerical calculations. 

In this study, we carry out the perturbation expansion 
to the fourth order in U for the self-energy S(k). The 
self-energy S(k) inu = U/t is expanded as 

E(k) = S {2) (k)u 2 + E^\k)u 3 + S (4 Hk)u 4 , (7) 

where we neglect the first-order term, that is, the Hartree 
term. This term provides only a constant shift, and can 
then be included by the chemical potential shift. All di- 
agrams that appear in the fourth-order perturbation are 
illustrated in Figs. 1 and 2. By numerically calculating 
the coefficient of each order term, S^(k), S^(k) and 
£( 4 \k), we investigate the convergence of the pertur- 
bation expansion in U. Each term is evaluated by the 
following set of equations. 

The bubble and ladder terms included in several dia- 
grams are given by 

Xo(q) = -££o(fc)So (£-<?), (8a) 
k 

Mq)= £So(fc)So(g-fc), (8b) 
k 

where q = (q,iv n ), and v n = 2nirT is the boson Mat- 
subara frequency. With these terms, the second- and the 
third-order terms are represented as 

£ {2) (k)= 5>o(«)0o(fc-«), (9a) 




id) (e) (/) 




Fig. 2. Diagrams of the fourth-order term E^(k) in self-energy, 
(a) — (c) correspond to £^ { (k), which is the reducible self-energy 
including Green's function with the second-order self-energy; 
(°0 — (/) -^rpaW' wn i cn i s the self-energy mediated by one 
boson fluctuation; (g) — (i) ^^(fc), which denotes the vertex 
correction for one boson fluctuation itself; (j) — (I) S^ 2 (k), 
which denotes the vertex correction with two boson fluctuation 
crossing. 

= -J2Mq)Go(q-k), (9b) 

and 

^ (3) (fc)=£{xo(<z) 2 So(fc-<z) 

On the other hand, twelve terms in the fourth-order per- 
turbation are classified into four groups. Each group con- 
tains three terms, which are exactly equivalent at half- 
filling n = 1, (See Appendix) 

S^(k) = S^+S^^ + S^l^ + S^k). (11) 

Here, E^ t (k) is the reducible self-energy including 
Green's function with the second-order self-energy. The 
corresponding diagrams are denoted by (a) — (c) in Fig. 2. 
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■^rpaCO * s ^ e self-energy mediated by one boson fluctu- 
ation, which is represented by (d) — (/) in Fig. 2. These 
diagrams, which are included in the FLEX, display fluc- 
tuations of the longitudinal mode, transverse mode and 
ladder diagram, respectively. E^ifo) is the vertex cor- 
rection for one boson fluctuation itself (type I), which is 
denoted by (g) — (i) in Fig. 2. These include, respectively, 
the longitudinal mode revised by itself, the ladder dia- 
gram revised by the transverse mode, and the transverse 
mode revised by the ladder diagram. S^l 2 (k) is the ver- 
tex correction with two boson fluctuations crossing (type 
II), which is denoted by (j) — (I) in Fig. 2. (j) denotes the 
longitudinal mode crossing with itself, and both (k) and 
(I) express the transverse mode crossing with the ladder 
diagram. Each term is represented as 



(12a) 



(12b) 



4tliW - £{-(xi(?) + Xi(«))0o(* - ?) 

+ MQ)Qo(Q-k)}, 

zltUk) = £{A(fc, q )xv(q)Gv(k q) 

+ 2A'(k,q)Mq)Go(q-k)}, 

where 



Gi(k) 


= Go(k)£ (2) (k)Go(k), 


(13) 


xi(q) 


= J2MP,q)Go(p)G ( P -q), 


(14a) 


Mq) 


p 

= J2^(p,q)Gv(p)Gv(q- P ), 


(14b) 


x'M 


p 

= J2^"(P,q)Go(-p)Go(q-p), 
p 


(14c) 


A(k,q) 


= £xo(fc'-*O0o(fc , )0o(fc , -g), 

k' 


(15a) 


A'(k,q) 


= J2xo(k' -k)g (k')go(q-k'), 

k' 


(15b) 


A"(k,q) 


= Y l Mk , -k)g {k , )g {k , -q). 


(15c) 



By computing all these terms, we examine the perturba- 
tion expansion to the fourth order. However, it is very 
hard to calculate these terms, particularly, the vertex 
correction terms. We need some technical procedures. 
Before we proceed to numerical results, let us introduce 
them. 

3. Numerical Recipes 

In this section, we introduce some techniques of nu- 
merical calculations in the perturbation expansion. For- 
tunately, most of the terms in the self-energy possess 
the convolution form. They can be easily evaluated with 



the use of the Fast Fourier Transform (FFT) algorithm. 
However, since <?o(fc, iu n ) is not periodic as a function of 
u n , applying the FFT to the Matsubara frequency sum 
yields extra numerical errors. Instead, we start numerical 
calculations by (/o(fe, r), which is defined by 

g (fc,r)--(l-/(a))e-« fcT (16a) 

w\e-e* T (&>o). v ; 

The final expression is convenient in numerical calcula- 
tions due to few errors. The Fourier transform is defined 
by 

G a (r,r) = ^J2 g o(k,r)c- ik - r . (17) 

k 

With this relation, the bubble and ladder terms are ex- 
pressed in the forms 

Xo(q,ifn) = ^2in4d,TXo(r,T)e iqr , (18a) 

r 

= mt$dTg (r, r)0 o (r, P - r)e iqr , (18b) 

r 

MQM =J2' mt o dT Mr,r)e iqr , (18c) 

r 

= J2^t P drGo(r,T) 2 e [ ' ir , (18d) 

r 

where 



\u(,\dT ■ ■ ■ ~ — 



Here, N T is the number of meshes along the imagi- 
nary time axis. In this case, the cutoff of Matsubara 
frequencies is ttTN t . In order to reduce numerical er- 
rors, we need to set this cutoff value large enough. Al- 
though we can easily transform r into q in the above 
equations, we require some care in the integral over r. 
Xo(q,t) and (j>o{q, r) abruptly decrease far from r = 
or (3. In addition, </>o(<7,t) has two different limiting val- 
ues at r = ±0 or {3 ± 0, and then does not strictly 
match with the FFT algorithm. In order to avoid this 
problem, wc use interpolation at around r = and 
(3. For a given q, we carry out fitting by the function 
/(t) = (ao+air) exp(-a 2 r) + (fe +fe 1 r') cxp(-fe 2 T'), with 
t 1 — (3 — t. Then, 5(f> a (q 7 t) = <fro(q, r) — /(r) is a smooth 
function at r = and (3, where its value and slope almost 
vanish. We can apply the FFT with high precision, / (t) 
itself can be easily integrated analytically. Because of this 
careful treatment, Xo{Qi' lu n) and 4>o(q, iv n ) recover suit- 
able v n dependences in the high-frequency region, and in 
the particle-hole symmetric case, satisfy the exact rela- 
tion Xo(q + Qi i^n) = <fio{q, w n) within numerical errors, 
where Q = (tt,7t). (See Appendix) In Fig. 3, we illus- 
trate the v n dependences of Xo(<Z> ll/ n) and $o{q, 'w n ) at 
q = (tt, tt/2) for (T, fi) = (O.li, -OAt) on the logarithmic 
scale as a sample. One can see that the real and imagi- 
nary parts of them have suitable behaviors in the high- 
frequency region (respectively, is^ and v\ dependences), 
owing to the interpolation at around t = and (3. This 
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where 



Fig. 3. Xo(q, i^n) and 4>o(q, \v n ) at q = (w,n/2) in the high- 
frequency region. cj>' Q and <f>Q represent the real and imaginary 
parts of </>o, respectively. A subscript (00) corresponds to the 
calculation with (without) the interpolation, xo an d <t>0 show 
suitable v n dependences in the high-frequency region, while xoo 
and </>oo are curved owing to errors caused by applying the FFT. 
This is the case of generic q. The numerical calculation is carried 
out for (T, fj.) = (O.lt, -OAt) with N x N x N T = 64 x 64 x 1024. 



is the case for generic q, although xo(q, at q = (0, 0) 
vanishes for v n ^ 0. 

We next proceed to the calculation of the self-energy. 
The second- and third-order terms have a convolution 
form and arc evaluated as 



(20a) 



-<Mr)So(r,/S-T)}e'", 

where 

X (r) = J2xo(q) 2 c- i ' ir , (21a) 
<Po(r)=J2Mq) 2 ^ lqr - (21b) 



Here, the sum over r = (r, r) denotes 

£ = £-^-#££- 

r r* r r r 

In the particle-hole symmetric case, the third-order term 
vanishes within numerical errors owing to the interpola- 
tion, since the first term in cq. (20b) cancels out the sec- 
ond term. Also in the fourth-order terms, the self-energy 
correction term eq. (12a) and the RPA term eq. (12b) 
possess the convolution form, and then are represented 
as 



4 4 if(fc) = £{2Xo(r)£i (r,r) 



S$A(*)=£{2*i(r)0o(r,r) 

+ *i(r)Qo(r,0-T) 



}e ifcr , 



(22) 



(23) 



^i(r)-^xo(g) 3 c- i( "\ 
<lM'-> = >X«) 3 e-^. 



£< 

9 



(24a) 
(24b) 



All these terms in the convolution form can be calculated 
with the FFT algorithm. However, the vertex correction 
terms E^ xl (k) and i7^ 2 (fc) cannot be written in the 
convolution form. Their evaluation requires very difficult 
computational effort. A typical form in the vertex correc- 
tion terms is that of X~l(q)i which included in -S^txi (k). 
We need an efficient way of computing it. Let us next 
introduce the technique adopted in this paper. 

First of all, by the Fourier transform of XoiPi —P2), wc 
decomposes Xi(s) m to the sum of the convolution form 



xi(q) = £ {xo(pi -P2) 

x Go(pi)Go(pi 



Pl,P2 



q)Go(p2)Go(p2 + <?)} 



= £ £{xo(r)e^-^ 

PUP2 " xGo(pi)Go(pi+q)Go(P2)Go(P2+q)} 

= £xoW{£e ipir ao(pi)ao(pi+9)} 

r Pl (25) 

x{Y,e- ip2r Go(P2)Go(P2+q)}. 

In the final expression, for a given r, the sum over p\ 
and P2 can be evaluated with the FFT algorithm. How- 
ever, the calculation time, which rapidly increases with a 
square of M = N 2 N T , does not change yet. It is almost 
impossible to do it within a practical machine time. For- 
tunately, we can restrict the sum over r to a number 
much less than N 2 , considering a characteristic property 
of Xo( r , t). We next demonstrate it. 

In order to obtain the proper perspective and reduce 
numerical errors, wc classify the sum over r with the use 
of the space group Ci v symmetry. xo{ r , T ) has the same 
value at any r with the same distance \r\. There are eight 
equivalent points in generic r. (four equivalent points on 
the x- and y-axis and the boundary, and so on.) The sum 
of eight exponential functions with the same distance \r\ 
are transformed into the sum of eight C\ v basis functions 
expressed as 



£< 



MP1-P2) 



£ 0n(Pl)0n(P 2 ), (26) 



|r|,n 



where eight C^ v basis functions (f> n (p) at generic r = 
are classified with irreducible representations (IR), 
as shown in Table. I. Precisely speaking, the subscript n 
in <p n {p) should be |r| and n, where n denotes a kind 
of IR. To simplify the expression, however, we do not 
explicitly write |r| as a subscript of functions, hereafter. 
Instead, to explicitly show the expansion in C± v basis 



6 J. Phys. Soc. Jpn. 



Table I. Ca v basis functions <p n {p) at r = 
n IR Basis functions 

1 Aj \/2(cos(ip x )cos(jpy) + cos(jp x )cos(ip y )) 

2 A2 V2(sm(ip x )sm(jp y ) — s'm(jp x )sin(ip y )) 
V2~(cos(ip x ) cos(jp y ) - cos(jp x ) cos(ip y )) 
V2(sm(ip x ) s'm(jpy) + s'm(jp x ) sin(ip y )) 

2sin(ip^) cos(jp y ) 
2s'm(ipy) cos(jp x ) 
2s'm(jp x ) cos(ip y ) 
2sm(jp y ) cos(ip x ) 
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E 



functions, we denote Xo{ r ,i~) by x« r - Thus, 
Xi(q) = E Xnr{£^(Pi)0o(pi)ao(pi+?)e iwiT } 

\r\nr _ Pi 



= £ X nT fnr{q)fn-r{q), 
\r\nr 

where 

fnr(q) = £ 4>n{Pl)Qo{pi)Qo{P! + <z)c^ 1+ ^ 



Pi 



J2^(p^(Pi+qy { ^ + - )T , 



(27) 



(28a) 



Gn(Pl) = <t>n(Pl)Go(pi)- 



(28b) 

Since /„ t (ij) is in the convolution form, it is evaluated 
by calculating the Fourier transform as 

fnr(q) =J2^( r ^ T ^o(ri,-r 1 -t)c- [ ^. (29) 



As will be shown later, we take at most 20 points with 
large contributions for the sum over \r\ in cq. (27). 
This can be carried out using recent supercomputers. 
In this way, we can compute the vertex correction terms 

Xvtxi CO and ^vtx2(^)- Finally, the type-I vertex correc- 
tion term is given by 



KZi (*) = £{-(*> W + x * W) Go ( r ' r ) e 



-$ 2 Wao(r,/3-r)c ifcr } 



(30) 



where 



X 2 (r)=J2^ e ~ iqr > ( 31a ) 

x 2 (r) = £xi(<z)e-^, (3ib) 

$ 2 (r)=E^(9K i9r - (31c) 
<? 

Xl(«) = £ X°r/nr(«)/n-r(«), (32a) 
\r\nr 

X'M = E €rfnr{q)fnr{q)\ (32b) 
| r | nr 

<M<0 - E T&T&Mfn-M, (32c) 

\r\nr 
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f'nM) = E Gn( P l)Go(q ~ Pl)e^~^ 

Pl (33) 

= E^( r i' r i)^( r i' T i- r ) ci9ri - 

Tl 

The type-II vertex correction term is given by 

lrlnTk '' q x X o(q)Qn(k - q)e-^ 

+ 2 X nr g n (k')go(q-k')e i ^> 

x Mq)Gn(q - k)e~ iuT } 

= E Y,Xnr^ T {fnr{q)Xn{q)Qn(k-q) 

lrlnT q +2fU q )Mq)Qn(q-k)} 

= EE^^e^ 1 (34) 

|r-|raT ri 

x {/ nr (ri)a„(ri,Ti)-2/; T (ri)g n (ri,/3-Ti)}. 
When <pn{p) is an even function (n = 1, 2, 3, 4), 

a„(r)=E^( fc K ifcr > (35a) 

k 

fnr(r)=J2fnr(q)X a (q)^ lqr , (35b) 

fUr) = J2fnr(q)Mq)c- iqr \ (35c) 

9 

and when 4> n (p) is an odd function (n = 5, 6, 7, 8), 

gn(r) = -iJ2Sn(k)e- ikr , (36a) 

/ nr (r)= iJ2fnr(q)Xo(q)e- i<!r , (36b) 

9 

tW^i^tf^of^. (36c) 

9 

Thus, we can calculate a set of equations for the self- 
energy in the fourth-order perturbation within a practi- 
cal machine time. 

Before we finish this section, let us show that the re- 
stricted r summation is efficient. In Fig. 4, we display 
Xo(f , t = 0) and (f>o{r,T = 0). It is clear that they have 
large magnitudes at around r = (0,0). We can approx- 
imately represent the r summation by several points at 
around r — (0,0). In this paper, we restrict the r sum- 
mation to several points around r = (0, 0) at which the 
magnitude of functions exceeds 10~ 4 . This corresponds 
to the fact that Xo(q) and 0o(<?) are approximated by 
the smeared functions x,o(q) and 0o(?), as illustrated in 
Fig. 5. One can see that the smeared functions restore the 
overall features. In this case, we can confirm the quality 
of the approximation by adopting the above-mentioned 
method for the second-order self-energy. 

£ {2) (k)=J2xo(k-k>)g Q (k') 

k' 

= E E X a nMk)Mk')e [uJT e-^' T go(k') 

k' \r\nr 

= Ex°^n(fe)e w E^( fc >" iW ' r 

|r|nr k' 
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Fig. 4. xo(f,-r = 0) and </> (r, r = 0) for (T, /Li) = (0.U, -0.4t). 
They have large magnitudes only at around r = (0, 0). 



\r\nr 

= E^( fe )E^ft(0,r)e iUT . (37) 

|r| r 

The last equality results from the fact that Q n (Q, r) is 
finite, only when <p n (k) has symmetry (i.e., n = 1). 
Figure 6 exhibits a comparison between the direct cal- 
culation of the convolution form E^ 2 '(k) and the ap- 

(2) 

proximate value by the restricted r summation Sa (k). 
We can see little difference in quantity between these 
two self-energies on this scale. The relative errors are 
\S^ 2 \k) - S^\k)\/\S^\k)\ < 6 x 1(T 3 . Thus, this ap- 
proximation proves to be efficient. This fact can also be 
verified in the third-order self-energy. 

4. Numerical Results 

Using the technique introduced in §3, we calculate all 
terms up to the fourth order in self-energy. In this paper, 
we carry out practical numerical calculations for M = 
N 2 N T = 64 x 64 x 1024. We set T = 0.1* except in 
the case of the discussion of temperature dependence. 
First, we investigate the convergence of the perturbation 
expansion. Next, we demonstrate the behavior of single- 
particle quantities, such as the DOS. 

4-.1 Convergence of perturbation expansion 

As mentioned in §1, in the perturbation expansion for 
the impurity Anderson model, physical quantities behave 
like the exponential function as a function of U. Thus, 
the coefficient of U n in the perturbation series rapidly 




Fig. 5. (Color online) xo(<J, 0) and <f>o{Q, 0) calculated using 
eq. (18), and Xoili 0) an d </>o(<?> 0) calculated using the restricted 
r summation for (T, /i) = (O.lt, — 0.44). Although the structure 
of the latter functions xo an d 4>o ls somewhat smeared, they well 
restore the overall features. The unit of the x- and j/-axes is it /a, 
where a is the lattice constant. 




COn/t 

Fig. 6. (Color online) Real and imaginary parts of the second- 
order term at k = (tt, 0) for (T, fi) = (O.lt, —OAt). We can see a 
small difference in quantity between the direct calculation of the 
convolution form in eq. (20a), (fe, iuin), and the approximate 
value by the restricted r summation, S a (k,iui n ) in eq. (37). 



decreases approximately in proportion to ~ 1/n!. Owing 
to this property, even if we truncate the perturbation 
expansion within a finite order, these physical quantities 
rapidly approach the exact values with increasing cutoff 
order. In the impurity Anderson model, the fourth-order 
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perturbation possesses a sufficient accuracy in a moder- 
ate correlation regime. Here, with these well-known facts 
in mind, let us examine the convergence of the perturba- 
tion expansion for the self-energy in the lattice system. 
We first study the half-filled case, and next, the doped 
case, and then, discuss the validity of the perturbation 
expansion truncated within a finite order. 

Half-filled case At half-filling n = 1 , the chem- 
ical potential fj, = 0. In this case, the system possesses 
particle-hole symmetry. The FS is just at the AF Bril- 
louin zone boundary. Owing to the complete nesting with 
Q = (tt, 7r), the spin susceptibility x(QjO) Wlu diverge 
at the Neel temperature Tn. The system undergoes the 
phase transition to the AF phase, although Tn — > by 
the Mermin-Wargner theorem in exact two-dimensional 
systems. The FL state becomes unstable in the ground 
state. For T > Tn, however, the system can stay in the 
FL state, if we consider Anderson's continuation prin- 
ciple. In this section, we explain that the perturbation 
expansion can be valid for T > O.lt at least, and then, 
for large U > 6i, the FL state seems to break down par- 
tially. The remarkable features become clearer by calcu- 
lating single-particle quantities, such as the DOS. 

In Fig. 7, we display u> n dependence of the fourth-order 
term (fc, iu) n ) at the two points on the FS, fc = (tt, 0) 
and (7r/2,7r/2), and at the Tma point far from the FS, 
fc = (0,0). The real part KeS^ (k, iu n ) is equal to zero 
on the FS, and its magnitude is the largest at the Fma 
point, while the imaginary part ImS^ (k,iuj n ) has the 
largest magnitude at k = (tt, 0). The behavior of the real 
part on the FS indicates that the FS is not deformed by 
the electron correlation, since the real part at uj n — » 
limit is connected with an energy shift after the analytic 
continuation on the real axis. In Fig. 8, we compare the 
largest value in the fourth-order term with that in the 
corresponding second-order term. The third-order term 
vanishes due to particle-hole symmetry. (See Appendix) 
The overall behavior of ImS^ (k, iuj n ) is almost inde- 
pendent of the wave number k. At a glance, one can 
find that the fourth-order term is much smaller than the 
second-order term. This means that the perturbation ex- 
pansion converges rapidly. In Fig. 9, we evaluate the ratio 
of the fourth-order term to the second-order term, quan- 
titatively. For simplicity, we display only the ratio of the 
imaginary part, because it is larger than the real part, 
and never crosses 0. In addition, the behavior is impor- 
tant for the behavior of the quasiparticle on the FS, since 
the imaginary part at lo„ — > limit is connected with the 
mass renormalization factor after the analytic continua- 
tion. In Fig. 9, the ratio is large on the FS, and becomes 
the largest at fe = (n, 0). Although the value increases 
with smaller uj n , the maximum value is ~ 0.07. If we ex- 
pect the same degree of convergence in the perturbation 
expansion as ~ 1/n! in the impurity case, then the ratio 
has to be of the order of ^/^ = 1/12 ~ 0.08. Thus, in 
this case, the convergence of the perturbation expansion 
is considered good. At lower temperatures, however, the 
maximum value becomes larger than 1/12 as shown in 
the inset of Fig. 9. Therefore, for T < O.lt, the pertur- 
bation expansion may break, or at least we need higher- 
order terms. In order to clarify this point, we need to 
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Fig. 7. (Color online) Real and imaginary parts of the fourth- 
order term Z< 4 ' (fc, iw n ) at fe = (0,0), (tt, 0) and (ir/2,7r/2) at 
T = O.lt in the half- filled case. The imaginary part at fc = (7r,0) 
in the low-frequency region provides the largest magnitude. 
It is increasing roughly according to the logarithmic function, 
0.00034 log(0.148o; n ), represented by the dotted line. 
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Fig. 8. (Color online) Comparison between the second-order term 
Z' 2 ) (fc, iuj n ) and the fourth-order term Z' 4 ) (fc, iu)„) in the half- 
filled case. We display the real parts at fc = (0, 0) and the 
imaginary parts at fc = (n, 0). At a glance, one can see that 
Z^ 4 ) (fc, iui n ) is much smaller than Z' 2 ' (fc, ko> n ). 



evaluate the sixth-order term. However, this is difficult 
to perform. On the other hand, in the high-frequency re- 
gion, the ratio is almost zero. The convergence of the per- 
turbation expansion is very good. Thus, for the present 
case, we can expect that the perturbation expansion is 
applicable for T > O.lt, since the coefficient of each order 
term becomes rapidly small as in the impurity case. 

Next, let us separately examine the contributions of 
each type of term in the fourth-order self-energy. In 
Fig. 10, we show the imaginary part of each term at 
k = (it,0) and (0,0). Only the RPA term has a large 
negative value, and the other terms behave to compen- 
sate for it. Although its magnitude at fc = (tt, 0) con- 
tinues to increase in the low-frequency region, it turns 
to decrease at generic fc points far from the FS. As a 
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Fig. 9. (Color online) Ratio of ImZ 1 * 4 ) (fc, to lmE^(k,iuJ n ) 
at fc = (0,0), (tt,0) and (tt/2, n/2) in the half-filled case. The 
ratio is the largest at fc = (it, 0) at any frequencies, and the 
maximum value is smaller than jr/^r = 1/12 ~ 0.08. In particu- 
lar, the ratio is negligible in the high-frequency region. The inset 
denotes the ratio at fc = (tt,0) for T = 0.02t, 0.05t and O.lt as 
a function of log(cj n /t). The largest value at w n = ttT quickly 
increases over 1/12 ~ 0.08. 
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Fig. 11. (Color online) Magnitude of the imaginary part of the 
self-energy up to the second-order, u 2 lmS^ (fc), and up to the 
fourth order, baE(k) = u 2 lm£( 2 )(k) + u 4 lmEW(k), at fc = 
(•7r,0) and (7r/2,7r/2) for U = 13t on the logarithmic scale. The 
unit of the vertical axis is u 2 t = U 2 /t. The self-energy to the 
fourth order is closer to the atomic limit behavior, U 2 /Aiuj n . 
We can expect that most of contributions from the higher-order 
terms arc confined in the shaded area in the low-frequency region. 
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Fig. 10. (Color online) Imaginary part of each type of term in 
the fourth-order self-energy at fc = (it, 0) and (0, 0) in the half- 
filled case. The large negative value of the RPA term is cancelled 
out by the other terms almost completely in the high-frequency 
region. 



whole, the large value in the low-frequency region be- 
comes one-order smaller due to the compensation. On 
the other hand, in the high-frequency region, the self- 
energy correction term and the type-II vertex correction 
term are suppressed abruptly. The RPA term and the 
type-I vertex correction term have a long tail with op- 
posite signs. These terms cancel out each other almost 
completely in the high-frequency region. Therefore, the 
fourth-order term as a whole vanishes rapidly in the high- 
frequency region. Thus, the behavior of self-energy in the 
high-frequency region is dominated by the second-order 
term. This indicates that the simple perturbation expan- 



sion up to the fourth order keeps the exact atomic limit 
behavior U 2 /4iui n in the high-frequency region. 30 This is 
different from the approximate method involving the re- 
summation of the specific diagrams, such as the FLEX. 
This fact is important in realizing the Mott-Hubbard 
character. The failure of the incoherent Hubbard peak 
structure in the FLEX comes from this fact partially. 

Let us consider the behavior in the low-frequency re- 
gion in details. The increase in the ratio at (ir, 0) orig- 
inates from the increase in the fourth-order term itself. 
It roughly behaves like log(w„) as shown in Fig. 7. Such 
behavior is different from the conventional FL behav- 
ior. Generally, the imaginary part of self-energy in the 
FL state has a tendency to decrease its magnitude in 
the low-frequency region such as that observed for the 
second-order term. This implies that the self-energy in- 
cludes a term proportional to iuj n in the expansion at 
around u) n — > 0. This is equivalent to existence of the 
weight of the FL quasiparticle. Thus, the singular behav- 
ior like log(a->„) in the fourth-order term is not the con- 
ventional FL behavior, and rather, is regarded as the pre- 
cursor into the Mott transition, although the singularity 
is weaker than l/iu> n in the Mott transition. In Fig. 11, 
we can clearly confirm the asymptotic behavior. It indi- 
cates the magnitude of the imaginary part of the self- 
energy up to the second order, u 2 lmS^ (k), and up to 
the fourth order, ImS(k) = u 2 lmS^(k) + u 4 ImS^(k), 
at fc = (tt,0) and (tt/2, tt/2) for large U = 13*. They 
asymptotically approach the atomic limit U 2 /4iuj n re- 
gardless of fc in the high-frequency region. For the large 
U = 13t, the fourth-order perturbation is closer to the 
atomic limit form than the second-order perturbation. 
This implies that the fourth-order perturbation can re- 
flect the Mott-Hubbard character more strongly, which 
has not been grasped in the second-order perturbation so 
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Fig. 12. (Color online) Fourth-order terms I7( 4 ) (fc, iui„) at fc = 
(0,0), (||7r, 0) and (±§7r, ±§tt) for (T, /i) = (0.U, -0.72t) in the 
doped case. The magnitude of the imaginary part is about 3 
times smaller than the maximum value at fc = (it, 0) in the half- 
filled case. It decreases in the low-frequency region, which is the 
conventional FL behavior. 
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Fig. 13. (Color online) Imaginary part of the nth order terms for 



n = 2, 3, 4 at fc = (^7r,0) in the doped case. One can see 
that the third- and fourth-order terms are much smaller than 
the second-order term. 



far. In fact, as indicated later, the DOS shows the pseu- 
dogap behavior like the Mott-Hubbard AF gap, and the 
spectral weight at fc — (tt, 0) does not indicate the quasi- 
particle peak at u) = any longer. Thus, we can consider 
that the fourth-order perturbation theory well describes 
the partial breakdown of the FL state. In order to more 
accurately describe the breakdown of the FL state, we 
require higher-order terms. We can easily imagine that 
such higher-order terms will become important in the 
shaded area in the low- frequency region in Fig. 11. The 
difference from the simple atomic limit is the remaining 
dispersive behavior in the low-frequency region. 

Doped case Next, let us consider a case except 

that for the half-filling. For instance, /i = — 0.72i cor- 
responds to n ~ 0.9 at U = 6t. The system possesses 
neither particle-hole symmetry nor the complete nest- 



Fig. 14. (Color online) Ratios of LmZ'f 3 ) (fc, iu) n ) and 
ImI7( 4 )(fc,iw n ) to ImEW(k,iuj n ) at fc = (0,0), (§§7T,0) 
and (^551", 351") in the doped case. The ratio of the third order 
is negative. This is because only the third-order term possesses 
positive values, as shown in Fig. 13. The ratios of the third and 
fourth orders are, respectively, smaller than 57 /jr = 1/3 and 
■h /ii = 1/12. The convergence of the perturbation expansion is 
considered very good. The inset denotes the ratios at (|^7r, 0) 
for T = 0.02t, 0.05t and 0.1* as a function of log(u> n /t). The 
maximum values do not become very large even for T < O.lt. 



ing. In this case, the FL state is considered more stable 
than that in the half-filled case. We can expect that the 
fourth-order perturbation is valid in a wider parameter 
region. In fact, we indicate that the perturbation expan- 
sion has better convergence, and thus, is applicable even 
at lower temperatures than that in the half-filled case. 

In Fig. 12, we display the oj n dependence of the fourth- 
order term at two points near the FS, fc — (§§tt, 0) and 
(j^TV, if 7r), and the Fma point. The magnitude of the 
imaginary part lmE^ 4 \k,iuj n ) is larger than that of the 
real part ReS^ (k,iuj n ), and has the largest value at 
k = (||7T, 0). The value is about 3 times smaller than 
the largest value at k = (it, 0) in the half-filled case. 
In addition, it turns to decrease in the low-frequency 
region. This is the conventional FL behavior, different 
from the behavior at k = (jr, 0) in the half-filled case. 
In Fig. 13, we illustrate the imaginary part of each or- 
der term at k = (||7T, 0). In this case, where the sys- 
tem is asymmetric with respect to the electron-hole, the 
third-order term (k) does not vanish and has a finite 
value. One can see that the third- and fourth-order terms 
are much smaller than the second-order term. This indi- 
cates that the perturbation expansion converges rapidly. 
In Fig. 14, we evaluate the ratios of ImS^ (k,iuj n ) and 
ImS^ (fc, iu! n ) to ImI7( 2 ) (fc, iu> n ) quantitatively. The dif- 
ference in sign between these ratios comes from the fact 
that only the third-order term possesses a positive sign, 
as shown in Fig. 13. This implies that the third-order 
term works to make the effective mass lighter. Thus, 
the third-order perturbation theory is insufficient for the 
effective mass. We need the fourth-order perturbation 
to obtain a large mass enhancement factor. In Fig. 14, 
we can see that the ratio of the fourth order is sup- 



J. Phys. Soc. Jpn. 



Full Paper 



Author Name 11 



pressed abruptly in the high-frequency region, and be- 
comes almost zero like that in the half-filled case. How- 
ever, strictly speaking, the ratio approaches a small but 
finite value, different from that in the half-filled case. 
This is because, in the doped case, the second-order term 
does not provide the exact form at the atomic limit, and 
generally, we require higher-order terms. 

Now, let us discuss the convergence of the perturba- 
tion series. If we expect the convergence of <~ 1/n! as dis- 
cussed in the half-filled case, the ratios of the third- and 
fourth-order terms have to be of the order of jr/^ = 1/3 
and 41/ 2i = 1/12, respectively. The actual ratios are 
smaller than those values. The convergence of the pertur- 
bation expansion is considered very good. In particular, 
the maximum ratio of the fourth-order term is smaller 
than that in the half-filled case, and does not become 
very large even at lower temperatures, as shown in the 
inset of Fig. 14. This indicates that the perturbation ex- 
pansion is valid even for T < O.lt, and the FL state 
is more stable than that in the half-filled case. This is 
because the doped system possesses neither the strong 
AF fluctuation by the complete nesting nor the Mott- 
Hubbard transition. 

However, we require some cautions in this consider- 
ation. First, we can take poU as the expansion param- 
eter, not the present U/t, where po is the DOS at the 
Fermi level. In this case, the ratio in the half-filled case 
becomes relatively smaller than that in the doped case, 
since po is largest in the half-filled case. Then the conver- 
gence of the perturbation expansion may almost be the 
same regardless of the electron filling. Second, we cannot 
simply apply the convergence criterion to the perturba- 
tion expansion in the present calculation. This is because 
the third-order term is relatively small in the vicinity of 
the half-filling since it completely vanishes just at half- 
filling. Therefore, we examined the ratio of the fourth- 
order term to the second-order term, not the ratio of con- 
secutive order terms. We need to evaluate higher-order 
terms to discuss the convergence of the perturbation sc- 
ries more accurately, although this is very difficult. 

Validity of the fourth-order perturbation From 

the above results, we would like to stress that the coef- 
ficient of the fourth-order term in the perturbation ex- 
pansion in U is very small as compared with that of the 
second-order term. For T > O.lt in the half-filled case, 
the maximum ratio is smaller than ~ 1/12 (= jf/^f), and 
in the doped case, it remains smaller even for T < O.lt. 
Such behavior of the coefficients is very similar to that 
of the specific heat coefficient in the impurity Ander- 
son model. 28 Also in that case, the odd-number order 
terms vanish in the symmetric case, and the coefficient 
of each nth-order term rapidly decreases roughly in pro- 
portion to 1/n!. The smallness of the fourth-order coef- 
ficient in the lattice system partially proves such good 
convergence. Thus, we can expect that the perturbation 
expansion in U maintains good convergence in a wide pa- 
rameter space of (n, T) also in the lattice system, except 
for T < O.lt in the half-filled case. In order to further 
clarify the behavior of ~ 1/n!, we need to investigate the 
coefficient of the sixth-order term, at least. It is one of im- 
portant future works, but difficult to perform. However, 



we consider that such good convergence in the pertur- 
bation expansion is inevitable, following the concept of 
the adiabatic continuation mentioned in §1. 29 As long as 
no phase transition occurs, the system connects adiabat- 
ically with the noninteracting system, and the physical 
quantities are analytic with respect to U . 

Now, let us consider the validity of the perturbation 
expansion truncated within a finite order. We examine 
the range of U where the second-order perturbation the- 
ory is valid, by comparing the second-order term and the 
fourth-order term, since the third-order term vanishes in 
the half-filled case. At T = O.lt in the half-filled case, 
the second-order perturbation is quantitatively valid for 
U < 3i ~ At from the ratio in Fig. 9. In the doped case, 
it holds even for T < O.lt from the inset in Fig. 14. For 
T > O.lt, the second-order perturbation is also valid for 
larger U. In the same way, to examine the validity of the 
fourth-order perturbation, we need the sixth-order term. 
In the present situation, we cannot discuss the validity 
quantitatively, since we do not estimate the contribu- 
tion of the sixth-order term. However, if the convergence 
of the perturbation expansion is good as we expected, 
and the coefficient of each nth-order term rapidly de- 
creases roughly in proportion to 1/n!, like that in the 
impurity case, then the fourth-order perturbation can be 
considered quantitatively valid for U < 5t <~ 6t from 
u 4 /4! ~ tt 6 /6!. From the smallness of the fourth-order 
coefficient and the concept of the adiabatic continuation 
mentioned above, we can expect that the expansion co- 
efficients behave like ~ 1/n! also in the lattice system, 
although we cannot guarantee it since we cannot esti- 
mate the correct contribution of neglected higher-order 
terms. In the following section, we display the DOS in the 
fourth-order perturbation for U < lOt at T = O.lt in the 
half-filled case. It shows a striking feature for U > 5t. We 
believe that the asymptotic behavior for large U in the 
fourth-order perturbation is well worth studying, even 
though it cannot be validated at present. Probably, the 
feature obtained in the fourth-order perturbation will be- 
come more distict by the inclusion of higher-order terms. 

4-2 Density of States 

We go on to a study of the DOS. In this section, 
we confine ourselves to the half-filled case. First, we re- 
view the DOS in the second-order perturbation and the 
FLEX. Then, we explain the result of the fourth-order 
perturbation. The DOS in each approximation is evalu- 
ated by carrying out k summation of Green's function 
with the corresponding self-energy, and then using ana- 
lytic continuation. That is, the DOS is given by 

p(w) = --lmg R (r = 0,u), (38) 

where Qr{t = 0,u) is the analytic continuation of 

Q{r = 0> n ) = ^2^T£(fc>„). (39) 
k 

The analytic continuation on the real axis is numerically 
calculated with the use of the Pade approximation. 31 

Second-order perturbation In Fig. 17(a), we ex- 
hibit the DOS in the second-order perturbation. The 



12 J. Phys. Soc. Jpn. 



Full Paper 



Author Name 



DOS (2nd) Large U limit 



DOS (2nd) 




Fig. 15. (Color online) DOS in the second-order perturbation at 
large u = U/t limit at T = O.li in the half- filled case, ft exhibits 
two peaks like the <5 function at around w ~ ±f/2 for U > 30t. 
This is consistent with the fact that the second-order perturba- 
tion is exact at the atomic limit t — > 0. In this case, the DOS is 
given by p(u) = \ (S(ui + + g(u, - %)) . 
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Fig. 16. (Color online) Second-order retarded self-energy at k = 
(it, 0) for U = 30t in the half-filled case. It asymptotically ap- 
proaches the atomic limit U 2 /iui in the high energy region. The 
real part of the self-energy crosses u) at around ui ~ ±U/2 = 
±15i, at which the imaginary part is negligible. Thus, the re- 
tarded Green's function possesses the real poles at the intersec- 
tions, and the DOS exhibits the S function peaks there. 



DOS (FLEX) 




DOS (4th) 




characteristic three-peak structure corresponds to the 
upper and lower Hubbard bands, and the coherent quasi- 
particle central peak. With increasing U, the structure 
becomes distinct. The width of the central peak corre- 
sponds to the quasiparticle bandwidth. With larger U, 
it becomes narrower, and the spectral weight is trans- 
ferred to the incoherent Hubbard peaks in the high en- 
ergy region. As shown in Fig. 15, at the large U limit, the 
quasiparticle weight almost vanishes, and the Hubbard 
peaks become similar to the 6 function. The DOS sub- 
stantially behaves like p{u) ~ \{5(uj + 2) + 8(u> — -j)). 
This can be understood from the behavior of the retarded 
self-energy Eji(k,u>) = u 2 (k,ui) for large U = 30t in 
Fig. 16. In this case, the real part crosses u> at around 



Fig. 17. (Color online) (a) DOSs in the second-order perturba- 
tion, (b) FLEX and (c) fourth-order perturbation at T = O.lt in 
the half-filled case. In (a), we can see the three-peak structure, 
which corresponds to the upper and lower Hubbard bands, and 
the coherent quasiparticle central peak. In (b), the pseudogap 
behavior at the Fermi level appears owing to the strong AF fluc- 
tuation, although the incoherent Hubbard structure is smeared 
and unclear. In (c), we find distinct incoherent peaks and the 
pseudogap structure for large U > 5i. This structure is consid- 
ered as the precursor of the Mott-Hubbard AF structure. 



oj ~ ±U/2 = ±15t, at which the imaginary part is neg- 
ligible. In other words, the retarded Green's function 
Qn(k,uj) possesses real poles at around uj ~ ±U/2. Such 
behavior in the second-order perturbation is reasonable, 
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not an artifact of the approximation, since it provides the 
exact form of self-energy at the large U limit, that is, the 
atomic limit (t — > 0). 30 This is a special circumstance in 
the half- filled case. In this case, in the atomic limit, the 
exact retarded self-energy is given by U 2 /4oj, and the 
Green's function possesses two poles at u> = ±Z7/2. In 
Fig. 15, we can see that the DOS almost restores such 
behavior for U > 30i, although the peak positions are 
still slightly larger. 

FLEX In Fig. 17(b), we display the DOS in the 

FLEX for various U values. The pseudogap occurs at the 
Fermi level u> — with a larger U owing to the strong AF 
fluctuation. On the other hand, the incoherent Hubbard 
bands in the high energy region are smeared and unclear. 
This is because the FLEX does not properly include the 
effect of the local correlation, as mentioned in §1. If we 
are reminded of our discussion of Fig. 11, this is related to 
the fact that the FLEX is inconsistent with the second- 
order perturbation in the high-frequency region, that is, 
it does not restore the atomic limit. 

Fourth-order perturbation Finally, the DOS in 

the fourth-order perturbation is illustrated in Fig. 17(c). 
For U < 3t, the behavior is almost the same as that in 
the second-order perturbation. It shows a qualitatively 
similar behavior even at U = At. For U > 5i, however, 
the shape is quite different from that in the second-order 
perturbation and the FLEX. It exhibits distinct incoher- 
ent peaks and the growth of the pseudogap at the Fermi 
level. The full-width of the DOS shrinks as compared 
with the second-order perturbation. This is considered 
as the band-narrowing effect caused by the correlation. 
In this case, the incoherent Hubbard bands are located 
at around the atomic limit positions already at U ~ W, 
where W = 8t is the bandwidth. This value is much 
smaller than U ~ 30t in the second-order perturbation. 
This is because the self-energy in the fourth-order per- 
turbation is closer to the atomic limit form, as shown 
in Fig. 11. The pseudogap develops with larger U, and 
the energy scale, ~ 2i, is much larger than that in the 
FLEX. Such a structure in the DOS is rather similar to 
that in the DCA, 7 although the peak structure of the 
incoherent part is sharper in the fourth-order perturba- 
tion. For U = W, the DCA has shown the formation 
of two coherent bands above and below to — owing 
to the strong AF correlation in addition to broad up- 
per and lower Hubbard bands. This characteristic Mott- 
Hubbard AF four band structure has been discussed by 
the quantum Monte Carlo method 7 ' 32 ' 33 and the exten- 
sion of DMFT. 11 Here, we have found that the DOS in 
the fourth-order perturbation also exhibits similar be- 
havior. 

Thus, we can consider that the fourth-order pertur- 
bation appropriately describes the asymptotic behavior 
into the Mott-Hubbard transition. In order to clarify 
such behavior, let us next investigate the behavior of 
self-energy at k points on the FS. 
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Fig. 18. (Color online) (a) Spectral weight, (b) the real part of the 
retarded self-energy and (c) the imaginary part for U = 6t at sev- 
eral k points on the FS shown in the inset of (a) in the half-filled 
case. The quasiparticle peak is well defined at fc = (7r/2,7r/2), 
far from which it gradually decreases, and then disappears at 
k = (it, 0). Correspondingly, the real part of the retarded self- 
energy at fc = (7T, 0) has a positive slope at ui = 0, at which the 
imaginary part indicates a large negative. Such behavior is quite 
different from the conventional FL behavior. 



4-3 Spectral weight and self-energy 

Figure 18 represents the spectral weight 

p(k,uj) = 



-lmQ R (k : uj) 



(40) 
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and the retarded self-energy £n(k,uj) at several fc points 
on the FS for U — 6t in the half-filled case. Note that, 
as shown in Fig. 18(a), the spectral weight at k — (tt, 0) 
does not exhibit the coherent quasiparticle peak struc- 
ture, although that at fe = (tt/2, tt/2) maintains the 
characteristic three peak structure. This indicates that 
the quasiparticle is not well defined at fe = (n, 0), al- 
though it still survives at k = (tt/2, 7r/2). In fact, the 
real part of S R (k,u) at k = (tt, 0) in Fig. 18(b) has a 
positive slope at u> = 0, at which the imaginary part 
in Fig. 18(c) has a large negative value. Such behav- 
ior is quite different from the conventional FL behavior 
observed at k = (tt/2, tt/2). Thus, the strong correla- 
tion breaks down the conventional FL quasiparticle at 
k = (tt,0). It may be rather better to consider that new 
two coherent bands are formed owing to the AF cor- 
relation, since the real part of the retarded self-energy 
has a negative slope at around u> ~ ±2i and the mag- 
nitude of the imaginary part is relatively small. In fact, 
the spectral weight seems to exhibit new hump struc- 
tures at around u> ~ ±t. 34 This is consistent with the 
formation of the Mott-Hubbard AF four-band structure 
mentioned above. 7 ' 11 ' 32 ' 33 In this case, however, far from 
k = (tt, 0), the structure is gradually suppressed, and the 
quasiparticle peak is restored. The pseudogap is large at 
k = (tt, 0), and vanishes at k = (tt/2, tt/2). This is differ- 
ent from the fully gapped behavior in the Mott-Hubbard 
AF structure. Although the DOS should be fully gapped 
at large U limit, in the fourth-order perturbation, even 
for larger U, the pseudogap is not open at k = (tt/2, tt/2). 
As discussed in Fig. 11, such gap behavior is related to 
the tendency toward an increase in the imaginary part of 
the self-energy ImS(k, iu>„) with a small uj n . ImS(k, iuj n ) 
at k = (tt/2, tt/2) does not show the tendency even for 
large U = 13t. The situation will be improved with the 
higher-order terms, since they are expected to possess a 
large contribution in the shaded area in Fig. 11. Thus, 
the pseudogap behavior in the fourth-order perturbation 
can be regarded as the precursor of the Mott-Hubbard 
AF four band structure. 

4-4 Mass enhancement factor 

The mass enhancement factor, which is the inverse of 
the mass renormalization factor z^, is given by the slope 
of the real part of the retarded self-energy at uj = 0, 

dRe£ R (k,oj) 



1 



dco 



(41) 



In Fig. 19, we illustrate the U dependence of Zk at 
k = (tt,0) and (tt/2, tt/2) in the half-filled case. In the 
second-order perturbation, we obtain z^ 1 = l + u 2 0.0166 
at fe = (vr,0) and z k 1 = 1 + w 2 0.026 at fe = (tt/2, tt/2) 
from the analytic continuation of (k, iu> n ). In the 
FLEX and the fourth-order perturbation, we carry out 
the Pade approximation for self-energy at each U. Con- 
cerning the quasiparticle at k — (tt/2, tt/2), which is al- 
ways well defined, z^ in the fourth-order perturbation is 
the smallest among the three approximations, and is very 
small for a large U. This indicates the formation of the 
quasiparticle with very heavy mass. The mass enhance- 
ment factor z^ 1 in the fourth-order perturbation is sev- 
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Fig. 19. (Color online) Mass renormalization factor Z& at fc = 
(•7r,0) and (7r/2,7r/2) in the half-filled case. In the second-order 
perturbation, zj, is given by 1/(1 + u 2 0.0166) at fc = (tt, 0) and 
1/(1 + u 2 0.026) at fc = (tt/2, tt/2). In the FLEX, z fe is more 
suppressed than that in the second order. In the fourth-order 
perturbation, z*. at fc = (it /2, it /2) has the smallest value. This 
means the formation of a quasiparticle with very heavy mass. 
On the other hand, z^, at fc = (it, 0) becomes negative at U > 
6t. This corresponds to the formation of the Mott-Hubbard AF 
structure. 



eral times larger than that in the FLEX. We can expect 
that it becomes larger at lower temperatures. Such be- 
havior is also obtained for the generic case except in the 
case of half-filling. Thus, the fourth-order perturbation 
theory can describe the quasiparticle with mass as heavy 
as that in heavy fermion systems. On the other hand, 
at fc = (tt,0), the fourth-order perturbation is qualita- 
tively different from the other two approximations. In the 
fourth-order perturbation, for U > 6t, z^ possesses a neg- 
ative value and the quasiparticle is not well defined. This 
corresponds to the formation of the Mott-Hubbard AF 
structure as mentioned above. Although we cannot deny 
that this may mean the breakdown of the perturbation 
expansion, it will be clarified by studying higher-order 
terms. 

4-5 Temperature dependence 

Finally, let us discuss the temperature dependence of 
the DOS in the fourth-order perturbation. In Fig. 20, we 
illustrate the DOS for various temperatures at U = 6t. 
Figures 20(a) and 20(b) show results in the half-filled 
case n = 1 and doped case n ~ 0.9, respectively. In 
Fig. 20(b), the quasiparticle peak at the Fermi level and 
the side peaks develop with decreasing temperature. The 
behavior of the former is consistent with the conventional 
FL state, although the latter are too enhanced due to 
the Pade approximation. Thus, in the doped case, the 
FL quasiparticles with heavy mass are formed with de- 
creasing temperature. On the other hand, in Fig. 20(a), 
the pseudogap develops at the Fermi level with decreas- 
ing temperatures. This corresponds to the formation of 
the Mott-Hubbard AF structure as discussed above. We 
cannot find such behavior for small U in the same tern- 
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Fig. 20. (Color online) DOSs at V = 6i (a) in the half-filled case 
n = 1 and (b) doped case n ~ 0.9. The origin of the vertical axis 
is shifted by 0.02 for each temperature. As shown in (a), with 
decreasing temperature, the pseudogap at the Fermi level devel- 
ops. This indicates the asymptotic behavior when the FL state 
breaks down. In contrast, in (b), we can see that the sharp peak 
at the Fermi level grows. Thus, in this case, a FL quasiparticle 
with heavy mass is formed. 



perature region. Thus, the fourth-order perturbation ex- 
pansion seems to describe the asymptotic behavior into 
the Mott-Hubbard transition. This unexpected result is 
very fascinating. This characteristic behavior will also be 
reflected in transport phenomena. It is very interesting 
to clarify whether the conductivity remains metallic. If 
it shows insulating behavior or a precursor, then we can 
study how to approach the metal-insulator transition. In 
order to clarify such behavior, we need further investiga- 
tion. All these will be investigated in our next study. 

5. Summary and Discussion 

In this study, we have investigated the Hubbard model 
on a two-dimensional square lattice by the perturbation 
expansion to the fourth order in the on-site Coulomb re- 
pulsion U. Numerically calculating all diagrams up to 
the fourth order in self-energy, we have examined the 
convergence of the perturbation series in the lattice sys- 
tem. The numerical techniques used here will be useful 
in calculating generic vertex correction terms, which can- 
not be reduced to their convolution form. In addition, 



we have reported the results of the DOS and the mass 
enhancement factor. Let us here summarize our results 
below. 

First, one of the most important results is the fact that 
the coefficient of the fourth-order term in self-energy is 
much smaller than that of the second-order term. The 
ratio is large at a small u) n , and the maximum values at 
T = O.Li are - 0.07 in the half-filled case and ~ 0.04 in 
the doped case. These values are consistent with the be- 
havior of ~ 1/n! for the nth-order term, which is almost 
the same as the behavior of the expansion coefficients in 
the impurity Anderson model. 28 Such good convergence 
in the perturbation expansion is considered as a general 
feature in the FL state from the concept of the adiabatic 
continuation. 29 The smallness of the fourth-order coeffi- 
cient partially proves such good convergence in the lattice 
system, although we cannot guarantee it since we can- 
not exactly evaluate neglected higher-order terms. Thus, 
we can expect that, as far as the FL state is stable, the 
lattice system also maintains good convergence, and the 
coefficients of higher-order terms are very small. In fact, 
the fourth-order term at a large u>„ almost vanishes. This 
is consistent with the fact that the FL state is fairly sta- 
ble at high temperatures. In addition, for T < O.lt, al- 
though the maximum ratio does not become very large 
in the doped case, it becomes large rapidly in the half- 
filled case. This is probably because the ground state at 
T = is the Mott-Hubbard AF state in the half-filled 
case, although the FL state is stable in the doped case. 
Thus, in the present lattice system, we can expect that 
the perturbation expansion in U keeps good convergence 
in wide parameter space of (n, T), except for T < O.lt in 
the half-filled case. 

Next, we evaluated the range of validities of the pertur- 
bation expansion truncated within a finite order. We can 
examine the range of U where the second-order perturba- 
tion theory is valid, by comparing the second-order term 
and the fourth-order term, since the third-order term 
vanishes in the half-filled case. From the above ratios, the 
second-order perturbation proves to be quantitatively 
valid for U < 3t ~ At. We cannot exactly determine the 
range of validity of the fourth-order perturbation, since 
we did not estimate the sixth-order term. However, fol- 
lowing the l/n!-like behavior mentioned above, we can 
expect that it is quantitatively valid for U < 5t ~ Qt 
from it 4 /4! ~ u 6 /6!. 

As for the behavior in the high-frequency region, the 
fourth-order term is abruptly suppressed owing to al- 
most perfect cancellation of contributions from each dia- 
gram. At half-filling, in particular, the second-order term 
provides the exact form of the self-energy in the high- 
frequency region. Thus, the self-energy to the fourth or- 
der keeps the atomic limit form [/ 2 /4kj„ in the high- 
frequency region. 30 This is important for the incoherent 
Hubbard peaks and the feature of the Mott transition. 
The failure in the FLEX comes partially from the fact 
that it is inconsistent with the atomic limit form in the 
high-frequency region. In the fourth-order perturbation, 
the self-energy for a large U more closely approaches the 
atomic limit form. This indicates that it can reflect the 
Mott-Hubbard characters more strongly, which have not 
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been grasped in the second-order perturbation so far. In 
fact, we found the asymptotic breakdown of the FL state 
in the strong correlation regime as the pseudogap behav- 
ior in the DOS. In order to clarify this feature more pre- 
cisely, we need to evaluate higher-order terms. We can 
expect that those terms decrease more rapidly in the 
high-frequency region, and most of the contributions are 
confined in the low-frequency region. 

Furthermore, we found interesting behaviors in the 
DOS. At half-filling, two features are noticeable. One 
is about the incoherent Hubbard structure in the high 
energy region. The peak structure is more remarkable, 
and the position is already very close to the atomic limit 
value co ~ ±U/2 for U ~ W, where W is the band- 
width 8t, although the second-order perturbation barely 
restores such behavior at around U ~ AW . Another 
is about the remarkable pseudogap phenomena at the 
Fermi level. With increasing U or decreasing temper- 
ature, the pseudogap develops. The energy scale <~ 2t 
is much larger than that of the simple AF gap ob- 
served in the FLEX. The pseudogap is rather consid- 
ered as the precursor of the formation of the Mott- 
Hubbard AF structure 7, n ' 32, 33 In this case, although the 
spectral weight at around k — (ir, 0) opens the pseu- 
dogap, the structure is gradually suppressed far from 
k = (it, 0), and the quasiparticle peak is restored at 
around k = (ir/2,ir/2). This is different from the con- 
ventional isotropic AF gap. Although the DOS should be 
fully gapped at a large U limit, the pseudogap remains 
close at k = (tt/2, n/2) in the fourth-order perturbation. 
This implies that higher-order terms are required at a 
larger U. Rather, such gapless behavior in the fourth- 
order perturbation seems to be consistent with the Fermi 
arc phenomena observed in the ARPES. 35 However, this 
should be discussed with a more realistic band struc- 
ture and electron filling. This is one of our future works. 
On the other hand, in the doped case, the DOS exhibits 
both a narrow quasiparticle peak at the Fermi level and 
an incoherent Hubbard structure. In this case, we ob- 
tain a very large mass enhancement factor. Thus, the 
fourth-order perturbation theory overall well explains the 
asymptotic behavior in the strong correlation regime. 

Finally, let us suggest several future works. The pseu- 
dogap behavior obtained in the fourth-order perturba- 
tion begins to develop in a rather high-temperature re- 
gion. In this case, it is very interesting to clarify whether 
the resistivity remains metallic. The study of transport 
phenomena based on the perturbation expansion is one of 
our future works. Another work is about AF and super- 
conducting transitions within the fourth-order perturba- 
tion. We can evaluate the superconducting (or AF) gap 
equation in the fourth-order perturbation theory. The 
preliminary results for the superconducting transition in- 
dicate that with increasing U, the eigenvalue of d-wave 
spin-singlet pairing increases abruptly in the weak cor- 
relation regime, and becomes optimal at around U ~ 5i, 
and then gradually decreases in the strong correlation 
regime. Such behavior is reasonable, and is also consis- 
tent with the result of recent variational Monte Carlo 
calculations. 36, 37 The increase of eigenvalue with increas- 
ing U originates from the increase of the attractive force, 
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and the decrease of eigenvalue for a large U limit results 
mainly from the increase of the mass enhancement fac- 
tor, namely, the renormalization effect. 21,22 Thus, the 
fourth-order perturbation theory throws light on investi- 
gations of the moderate correlation regime, which have 
been very difficult to perform so far. The phase diagram 
and physical properties will be discussed in our forth- 
coming paper. 
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Appendix 

In the symmetric case, the third-order term vanishes, 
and three terms contained in each type of the fourth- 
order term becomes equivalent. These exact relations can 
also be a test to check numerical calculations. Let us here 
prove these relations. 

The dispersion relation given in eq. (2) satisfies ek+Q = 
— £fe with Q — (7r,7r). At half-filling n = 1, the chemical 
potential fj, = 0. In this case, the single-particle Green's 
function possesses the property 

Go(k + Q) = G (k + Q,iuj n ) 
1 

iuj n + e k 
= -Qo(-k), 

where Q = (Q,0). The ladder diagram is related to the 
bubble diagram, 

4>(q + Q) = ^a (fc)g (g-fc + Q) 

k 

= -]TSo(fc)eo(fc-<?) 

k 

= x(q)- 

Thus, 

]T Mq) 2 Go(q - fc) = £ Mq + QfGoig -k + Q) 

q <i 

= - ^2xa(q)Qa(k - q), 

q 

and then, the third-order term is proved to vanish. 

Next, we consider each term of the fourth order. First 
of all, since the second-order self-energy is calculated as 

s^(k + Q) = -J2 Mq)Go(q -k-Q) 

q 

= -Y,Mq + Q)G°(q-k) 

q 

= -£xo(9)0o(«-fc) 
q 

= -S^\-k), 
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then we obtain 

Gi{k + Q) = -Gi(-k). 
Thus, in the self-energy correction term I^j(fc), 

- Mq)0i{q - fc) = -£>(<? + Q)Gi(q + Q-k) 
q q 

q 

and then, (a) — (c) in Fig. 2 have equivalent contributions. 
In the RPA term £$ A (k), 

- E ^o(Q)Oo (q - *0 = " I>o(« + Q)$o(q + Q-k) 

q q 

= J2x 3 o(q)Go(k-q), 
q 

and then, (d) — (/) in Fig. 2 have equivalent contribu- 
tions. Next, from A"(k + Q,q) = A(k, q), 

X 'M = A"(p + Q, q)Qo(-p - Q)Go(q -p-Q) 
= A(p,q)g (p)g (p-q) 

= xi(q), 

and from A'(k, q + Q) = — A(k, q), 

<t>i{q + Q) = A'(j>,q + Q)g {p)g (q-p + Q) 
= A(p,q)g (p)g (p-q) 

= xi(q)- 

Thus, in the type-I vertex correction term E^l^k), 

- E xi(q)Go(k -q) = -J2 x[(q)Go(k - q) 

q q 

= J2Mq)Go(q-k), 

q 

that is, (g) — (i) in Fig. 2 have equivalent contributions. 
Finally, in the typc-II vertex correction term S^ 2 (k), 

J2^'(k,q)Mq)Go(q-k) 

q 

= E A'(fc, q + Q)Mq + Q)G°(q -k + Q) 

q 

= J2Hk,q)Xo(q)Go(k-q). 
q 

Namely (j) — (I) in Fig. 2 have equivalent contributions. 
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